##  Project: CueAnon
##  Authors: Benjamin S. Noble and Taylor N. Carlson
## Code to replicate tables and figures for vignette experiment
## Figures 1-2, Appendix B.

library(margins)
library(RecordLinkage)
library(utf8)
library(texreg)
library(xtable)
library(rstatix)
library(kableExtra)
library(marginaleffects)
library(modelsummary)
library(ggpubr)
library(tidyverse)
options("modelsummary_format_numeric_latex" = "mathmode")

source('ggtheme_baselike.R')

# data from vignette experiment
wc <- read_csv('tass_small.csv') 

# --------------------------------------------------------------------------- #

# summary statistics
# wave
wc %>% 
  group_by(wave) %>% 
  summarise(n = n())

# proportion don't know qanon by wave 
wc %>% 
  group_by(wave) %>% 
  summarise(dk = sum(dk_q), n = n()) %>% 
  mutate(prop_dk = dk/n)

# proportion support qanon among those who know qanon
wc %>% 
  filter(dk_q == 0) %>% 
  summarise(sq = sum(support_q), n = n()) %>% 
  mutate(prop_s = sq/n)

# trust proportions
wc %>% 
  group_by(trust) %>% 
  summarise(n = n()) %>% 
  bind_cols(tot = nrow(wc %>% filter(!is.na(trust)))) %>% 
  mutate(prop_trust = n/tot)

# trust by party
wcpid <- wc %>% 
  mutate(rep = if_else(dem == 0 & ind == 0, 1, 0),
    trust_dummy = if_else(trust <= 2, 1, 0)) %>% 
  group_by(rep, trust_dummy) %>% 
  filter(!is.na(trust)) %>% 
  summarise(n = n())

wcpid %>% filter(rep == 1 & trust_dummy == 1) %>% pull(n)/wcpid %>% filter(rep == 1) %>% pull(n) %>% sum()
wcpid %>% filter(rep == 0 & trust_dummy == 1) %>% pull(n)/wcpid %>% filter(rep == 0) %>% pull(n) %>% sum()

# --------------------------------------------------------------------------- #

# Table B1, experiment balance
exp_bal <- wc %>% 
  group_by(exp_cond) %>% 
  select(exp_cond, weight, newsimport, pid7, ideology, female, age4, educ5, income, white, black, latin, asian, metro_area) %>% 
  summarise(across(newsimport:metro_area,  ~weighted.mean(.x, w = weight, na.rm = T))) %>% 
  t() %>% 
  as.data.frame() %>%  
  round(., 2) %>% 
  rownames_to_column()

## F-test on each covariate between all 3 groups
## significance indicates a difference among groups
comps <- matrix(NA, nrow = 12, ncol = 2)
i = 1
t <- summary(aov(newsimport ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1] # f-stat
comps[i,2] <- unlist(t[[1]][5])[1] # p-value
i = i + 1

t <- summary(aov(pid7 ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(ideology ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(female ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(age4 ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(educ5 ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(income ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(white ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(black ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(latin ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(asian ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

t <- summary(aov(metro_area ~ exp_cond, wc, weights = weight))
comps[i,1] <- unlist(t[[1]][4])[1]
comps[i,2] <- unlist(t[[1]][5])[1]
i = i + 1

# append f-stat and p-value to balance table
btab <- bind_cols(exp_bal[-1,], 
  comps[,1] %>% round(., 2), 
  comps[,2] %>% round(., 2)) 
# rename cols
names(btab) <- c('Variable', 'Negative', 'Neutral', 'Conspiracy', 'Ftest', 'pv') 
# reformat table
btab <- btab %>% select(Variable, Neutral, Negative, Conspiracy, Ftest, pv) 
 # add variable names for latex
btab$Variable <- c('News Importance (5 point)', 'Party ID (7 Point)', 
  'Ideology (5 Point)', 'Female', 'Age (4 Point)', 'Education (5 Point)', 
  'Income (20 Point)', 'White', 'Black', 'Latin', 'Asian', 'Metro')
btab <- btab %>% rename(`F-Stat` = Ftest, `p-Value` = pv)
# Table B1
kable(btab, 'latex', escape = FALSE, booktabs = TRUE)

# --------------------------------------------------------------------------- #

# Figure 1
trust_me_plots <- function(model, treat_model, var_neg, var_q, v_range, x_lab){ # input the treat - neutral model, qanon - negative model,
  #  the name of the negative variable, name of q variable, range of moderator (trust), name of the x label (trust)
  neg_margin <- margins(model, variables = var_neg, at = list(trust = v_range)) # get marginal effects for negative-neutral
  neg_margin2 <- summary(neg_margin, level = 1-(.05/12)) # bonferroni correction for 12 tests—4 levels of trust x 3 comparisons
  neg_margin2 <- as_tibble(neg_margin2) %>% mutate(comp = 'Neg - Neu') # save the effects and name comparison condition
  
  q_margin <- margins(model, variables = var_q, at = list(trust = v_range)) # get marginal effect for q-neutral
  q_margin2 <- summary(q_margin, level = 1-(.05/12)) # bonferroni correction for 12 tests—4 levels of trust x 3 comparisons
  q_margin2 <- as_tibble(q_margin2) %>% mutate(comp = 'Q - Neu')  # save the effects and name comparison condition
  
  treat_margin <- margins(treat_model, variables = var_q, at = list(trust = v_range)) # get AME and CI for treatment comparison
  treat_margin2 <- summary(treat_margin, level = 1-(.05/12)) # bonferroni correction for 12 tests—4 levels of trust x 3 comparisons
  treat_margin2 <- as_tibble(treat_margin2) %>% mutate(comp = 'Q - Neg') # save the effects and name comparison condition
  
  fulldf <- bind_rows(neg_margin2, q_margin2, treat_margin2) # put all effects into a single df
  fulldf$comp <- factor(fulldf$comp, levels = c("Neg - Neu", "Q - Neu", "Q - Neg")) # re-level factors into correct order
  
  trust_all <- ggplot(fulldf, aes(x = as.factor(trust), y = AME)) + 
    geom_point(aes(shape = comp, color = comp), position = position_dodge(w = .3)) + 
    geom_errorbar(aes(ymin = lower, ymax = upper, color = comp), width = 0, position = position_dodge(w = .3)) + 
    geom_hline(yintercept = 0, linetype = 'dashed') +
    labs(x = x_lab, y = 'Marginal Effect of Treatment') + 
    scale_color_grey(name = 'Comparison', start = 0, end = .6) +
    scale_x_discrete(breaks = c(1:4), labels=c('None \nAt All', 'Not Very \nMuch', 'A Fair \nAmount', 'A Great \nDeal')) +
    scale_shape_manual(name = "Comparison", values = c(15, 16, 17)) + 
    theme_baselike() + 
    coord_cartesian(ylim = c(-50, 10))
  
  return(trust_all)
}

ols_trust <- lm(therm ~ negative * trust + qanon * trust + wave, wc) # therm ~ treat-neutral x trust
ols_trust_treat <- lm(therm ~ qanon * trust + wave, wc %>% filter(neutral != 1)) # therm ~ qanon-negative x trust

# values of coefficients 
margins(ols_trust, variables = 'qanon', at = list(trust = 1:4))
margins(ols_trust, variables = 'negative', at = list(trust = 1:4))
margins(ols_trust_treat, variables = 'qanon', at = list(trust = 1:4))

me_trust_plot <- trust_me_plots(ols_trust, ols_trust_treat, 'negative', 'qanon', 1:4, '\nTrust in Media') # create me plots
# Figure 1
me_trust_plot 

# bootstrap intervals for high/low differences
set.seed(2022)
me_diff1 <- me_diff2 <- c()
for(i in 1:1000){
  wc_b <- wc %>% sample_n(size = nrow(wc), replace = T)
  ols_trust_b <- lm(therm ~ negative * trust + qanon * trust, wc_b)
  me1 <- margins(ols_trust_b, variables = 'negative', at = list(trust = c(1, 4)))
  me_diff1 <- c(me_diff1, summary(me1)$AME[2] - summary(me1)$AME[1])
  
  me2 <- margins(ols_trust_b, variables = 'qanon', at = list(trust = c(1, 4)))
  me_diff2 <- c(me_diff2, summary(me2)$AME[2] - summary(me2)$AME[1])
}
# bonferroni quantiles
quantile(me_diff1, c((.05/12)/2, 0.5, (1-(.05/12)/2))) # difference in me of negative for high/low trust
quantile(me_diff2, c((.05/12)/2, 0.5, (1-(.05/12)/2))) # difference in me of qanon for high/low trust

# without wave indicators
ols_trust_nw <- lm(therm ~ negative * trust + qanon * trust, wc) 
ols_trust_treat_nw <- lm(therm ~ qanon * trust, wc %>% filter(neutral != 1))
me_trust_plot_nw <- trust_me_plots(ols_trust_nw, ols_trust_treat_nw, 'negative', 'qanon', 1:4, '\nTrust in Media') 
# Figure B8 (top)
me_trust_plot_nw

# --------------------------------------------------------------------------- #

# Table B2
ols_trust_control <- lm(therm ~ negative * trust + qanon * trust + newsimport + 
  pid7 + ideology + female + age4 + educ5 + income + 
  white + black + latin + asian + wave, wc, weights = weight)
ols_trust_wave1 <- lm(therm ~ negative * trust + qanon * trust, wc %>% filter(wave==1))
ols_trust_wave2 <- lm(therm ~ negative * trust + qanon * trust, wc %>% filter(wave==2))

c_map <- c(
  'trust' = 'Trust in Media',
  'negative' = 'Negative',
  'qanon' = 'QAnon',
  'negative:trust' = 'Negative x Trust', 
  'trust:qanon' = 'QAnon x Trust',
  'wave' = 'Wave',
  'newsimport' = 'News Importance',
  'pid7' = 'Party ID (7)',
  'ideology' = 'Ideology (7)',
  'female' = 'Female',
  'age4' = 'Age (4)',
  'educ5' = 'Education (5)',
  'income' = 'Income (18)',
  'white' = 'White',
  'black' = 'Black',
  'latin' = 'Latin',
  'asian' = 'Asian',
  '(Intercept)' = 'Baseline'
)

modelsummary(list(ols_trust, ols_trust_control, ols_trust_wave1, ols_trust_wave2),
    coef_map = c_map,
    stars = TRUE,
    gof_omit = 'AIC|BIC|FE|Std.',
    output = 'latex',
    fmt = fmt_decimal(digits = 2))

# --------------------------------------------------------------------------- #

# Figure B1 and Table B3
# see above function as this works almost exactly the same
pid_me_plots <- function(model, treat_model, var_neg, var_q, v_range, x_lab){ 
  neg_margin <- margins(model, variables = var_neg, at = list(pid7 = v_range)) 
  neg_margin2 <- summary(neg_margin, level = 1-(.05/18)) # adjust for 3 treatments x 6 comparisons (i.e., we don't care about independents)
  neg_margin2 <- as_tibble(neg_margin2) %>% mutate(comp = 'Neg - Neu')
  
  q_margin <- margins(model, variables = var_q, at = list(pid7 = v_range)) 
  q_margin2 <- summary(q_margin, level = 1-(.05/18)) 
  q_margin2 <- as_tibble(q_margin2) %>% mutate(comp = 'Q - Neu')
  
  treat_margin <- margins(treat_model, variables = var_q, at = list(pid7 = v_range))
  treat_margin2 <- summary(treat_margin, level = 1-(.05/18))
  treat_margin2 <- as_tibble(treat_margin2) %>% mutate(comp = 'Q - Neg')
  
  fulldf <- bind_rows(neg_margin2, q_margin2, treat_margin2)
  fulldf$comp <- factor(fulldf$comp, levels = c("Neg - Neu", "Q - Neu", "Q - Neg"))
  
  pid_all <- ggplot(fulldf, aes(x = as.factor(pid7), y = AME)) + 
    geom_point(aes(shape = comp, color = comp), position = position_dodge(w = .3)) + 
    geom_errorbar(aes(ymin = lower, ymax = upper, color = comp), width = 0, position = position_dodge(w = .3)) + 
    geom_hline(yintercept = 0, linetype = 'dashed') +
    labs(x = x_lab, y = 'Marginal Effect of Treatment') + 
    scale_color_grey(name = 'Condition', start = 0, end = .6) +
    scale_shape_manual(name = "Condition", values = c(15, 16, 17)) +
    scale_x_discrete(breaks = c(1:7), labels=c('Strong \nDemocrat', 'Democrat', 'Lean \nDemocrat', 'Independent', 'Lean \nRepublican', 'Republican', 'Strong \nRepublican')) +
    theme(axis.text.x = element_text(angle=90, hjust=1)) +
    theme_baselike() + 
    coord_cartesian(ylim = c(-45, 10))
  
  return(pid_all)
}

# create ME plot for party id
pid_ols <- lm(therm ~ negative * pid7 + qanon * pid7 + wave, data=wc)
pid_ols_treat <- lm(therm ~ qanon * pid7 + wave, data = wc %>% filter(neutral != 1))
me_pid_plot <- pid_me_plots(pid_ols, pid_ols_treat, 'negative', 'qanon', 1:7, '\nParty ID') # plots
me_pid_plot # Figure B1

# without wave indicators
pid_ols_nw <- lm(therm ~ negative * pid7 + qanon * pid7, data=wc)
pid_ols_treat_nw <- lm(therm ~ qanon * pid7, data = wc %>% filter(neutral != 1))
me_pid_plot_nw <- pid_me_plots(pid_ols_nw, pid_ols_treat_nw, 'negative', 'qanon', 1:7, '\nParty ID') # plots
me_pid_plot_nw # Figure B8 (bottom)

# figure B8
nw_plots <- ggarrange(me_trust_plot_nw, me_pid_plot_nw, ncol = 1)

# Table B3
ols_pid7_control <- lm(therm ~ negative * pid7 + qanon * pid7 + trust + 
  newsimport + ideology + female + age4 + educ5 + income + white +
  black + latin + asian + wave, wc, weights = weight)
ols_pid7_wave1 <- lm(therm ~ negative * pid7 + qanon * pid7, wc %>% filter(wave==1))
ols_pid7_wave2 <- lm(therm ~ negative * pid7 + qanon * pid7, wc %>% filter(wave==2))

c_map2 <- c(
  'pid7' = 'Party ID (7)',
  'negative' = 'Negative',
  'qanon' = 'QAnon',
  'negative:pid7' = 'Negative x Party ID', 
  'pid7:qanon' = 'QAnon x Party ID',
  'wave' = 'Wave',
  'newsimport' = 'News Importance',
  'trust' = 'Trust in Media',
  'ideology' = 'Ideology (7)',
  'female' = 'Female',
  'age4' = 'Age (4)',
  'educ5' = 'Education (5)',
  'income' = 'Income (18)',
  'white' = 'White',
  'black' = 'Black',
  'latin' = 'Latin',
  'asian' = 'Asian',
  '(Intercept)' = 'Baseline'
)

modelsummary(list(pid_ols, ols_pid7_control, ols_pid7_wave1, ols_pid7_wave2),
    coef_map = c_map2,
    stars = TRUE,
    gof_omit = 'AIC|BIC|FE|Std.',
    output = 'latex',
    fmt = fmt_decimal(digits = 2))

# --------------------------------------------------------------------------- #

# Table B4
pattern_w1 <- c('smith', 'john')
pattern_w2 <- c('cunningham', 'john')

wc <- wc %>% 
  mutate(namerec = tolower(utf8_encode(cand_name))) %>% # encode and lowercase all names
  rowwise() %>% 
  # if wave == 1, get jw score for john smith, else for john cunningham
  mutate(jw_score = if_else(wave == 1, jarowinkler(namerec, 'john smith'), 
                            jarowinkler(namerec, 'john cunningham')), 
         # robustness, if wave == 1 & john or smith is mentioned ~ 1, 
         # else 0; same for w2 and john or cunningham
         my_score = if_else((grepl(paste(pattern_w1, collapse = '|'), namerec) & 
                               wave == 1) | (grepl(paste(pattern_w2, collapse = '|'), 
                                                   namerec) & wave == 2), 1, 0),
          my_score = if_else(is.na(namerec), NA_real_, my_score))

jw_ols_nowave <- lm(jw_score ~ negative + qanon, wc)
jw_ols <- lm(jw_score ~ negative + qanon + wave, wc) # jw score
jw_ols_control <- lm(jw_score ~ negative + qanon + wave + newsimport + ideology + female + age4 + educ5 + income + white +
  black + latin + asian + wave, weights = weight, wc) # jw score
js_logit <- glm(my_score ~ negative + qanon + wave, wc, family = 'binomial') # binary name recognition

name_coef <- c(
  'negative' = 'Negative',
  'qanon' = 'QAnon',
  'wave' = 'Wave',
  'newsimport' = 'News Importance',
  'trust' = 'Trust in Media',
  'ideology' = 'Ideology (7)',
  'female' = 'Female',
  'age4' = 'Age (4)',
  'educ5' = 'Education (5)',
  'income' = 'Income (18)',
  'white' = 'White',
  'black' = 'Black',
  'latin' = 'Latin',
  'asian' = 'Asian',
  '(Intercept)' = 'Baseline'
)

modelsummary(models = list(jw_ols_nowave, jw_ols, jw_ols_control, js_logit),
    coef_map = name_coef,
    stars = TRUE,
    gof_omit = 'AIC|BIC|FE|Std.',
    output = 'latex',
    fmt = fmt_decimal(digits = 2))

# --------------------------------------------------------------------------- #
# Figure 2 and Table B5 
wc_i <- wc %>% mutate(gop = if_else(pid7 > 4, 1, 0), low_trust = if_else(trust < 3, 1, 0))

# all models needed---includes full data with treatments and treatment comparison
# then split by pid and by trust
ideology_full <- lm(perc_ideo ~ negative + qanon, wc_i)
ideology_full_treat <- lm(perc_ideo ~ qanon, wc_i %>% filter(neutral == 0))
ideology_pid_int <- lm(perc_ideo ~ negative * gop + qanon * gop, wc_i)
ideology_pid_treat <- lm(perc_ideo ~ qanon * gop, wc_i %>% filter(neutral == 0))
ideology_trust_int <- lm(perc_ideo ~ negative * low_trust + qanon * low_trust, wc_i)
ideology_trust_treat <- lm(perc_ideo ~ qanon * low_trust, wc_i %>% filter(neutral == 0))

perc_ideo_mes <- bind_rows(
  # marginal effects for full sample
  bind_rows(
    avg_slopes(ideology_full, variables = c('negative', 'qanon'), conf_level = 1 - (.05/15)),
    avg_slopes(ideology_full_treat, variables = c('qanon'), conf_level = 1 - (.05/15))
  ) %>% mutate(term = c('Neg — Neu', 'Q — Neu', 'Q — Neg')),
  # marginal effects for party id split
  bind_rows(
    slopes(ideology_pid_int, variables = c('negative', 'qanon'), newdata = datagrid(gop = c(0,1)), conf_level = 1 - (.05/15)),
    slopes(ideology_pid_treat, variables = c('qanon'), newdata = datagrid(gop = c(0,1)), conf_level = 1 - (.05/15))
  ) %>% mutate(term = c(rep('Neg — Neu', 2), rep('Q — Neu', 2), rep('Q — Neg', 2))),
  # marginal effects for trust split 
  bind_rows(
    slopes(ideology_trust_int, variables = c('negative', 'qanon'), newdata = datagrid(low_trust = c(0,1)), conf_level = 1 - (.05/15)),
    slopes(ideology_trust_treat, variables = c('qanon'), newdata = datagrid(low_trust = c(0,1)), conf_level = 1 - (.05/15))
  ) %>% mutate(term = c(rep('Neg — Neu', 2), rep('Q — Neu', 2), rep('Q — Neg', 2)))
) %>% 
  mutate(mod_lev = c(rep('All',3), rep(c('Dem/Ind', 'Rep'), 3), rep(c('High Trust', 'Low Trust'), 3)),
    mod_lev = factor(mod_lev, levels = c('All', 'Dem/Ind', 'Rep', 'High Trust', 'Low Trust')),
    term = factor(term, levels = c('Neg — Neu', 'Q — Neu', 'Q — Neg')))

# figure 2
perc_ideo_plot <- ggplot(perc_ideo_mes, aes(x = mod_lev, y = estimate, color = term, shape = term)) + 
  geom_point(position = position_dodge(.3)) + 
  geom_errorbar(aes(x = mod_lev, ymin = conf.low, ymax = conf.high), width = 0, position = position_dodge(.3)) + 
  theme_baselike() + 
  geom_hline(yintercept = 0, linetype = 'dashed') + 
  labs(x = '', y = 'Less Conservative --- More Conservative', color = '', shape = '') +
  theme(legend.position = 'bottom', legend.direction = 'horizontal', panel.border = element_rect(colour = "black", fill=NA, linewidth=1)) + 
  geom_vline(xintercept = c(1.5, 3.5), linetype = 'solid') + 
  scale_y_continuous(breaks = seq(-1, 2, by = 0.5)) +
  scale_color_grey(name = 'Comparison', start = 0, end = .6) +
  scale_shape_manual(name = "Comparison", values = c(15, 16, 17))

ideo_c_map <- c(
  'negative' = 'Negative',
  'qanon' = 'QAnon',
  'low_trust' = 'Low Trust in Media',
  'gop' = 'Republican',
  'negative:gop' = 'Negative x Republican',
  'negative:low_trust' = 'Negative x Low Trust in Media',
  'gop:qanon' = 'QAnon x Republican',
  'low_trust:qanon' = 'QAnon x Low Trust in Media',
  '(Intercept)' = 'Baseline'
)
# Table B5
modelsummary(models = list(ideology_full, ideology_pid_int, ideology_trust_int),    
    coef_map = ideo_c_map,
    stars = TRUE,
    gof_omit = 'AIC|BIC|FE|Std.',
    output = 'latex',
    fmt = fmt_decimal(digits = 2))

# --------------------------------------------------------------------------- #
# robustness tests
# additional mturk testing
mturk <- read_csv('mturk_followup_small.csv')

# Figure B2
win_quote_trust <- lm(fav ~ win_quote_negative * trust + win_quote_qanon * trust, mturk) # favorability quote, trust
win_quote_trust_treat <- lm(fav ~ win_quote_qanon * trust, mturk %>% filter(win_quote_neutral != 1)) # favorability quote, trust treatment comparison
win_quote_trust_plot <- trust_me_plots(win_quote_trust, win_quote_trust_treat, 'win_quote_negative', 'win_quote_qanon', 1:4, 'Trust in Media') + 
    coord_cartesian(ylim = c(-65, 10)) # plots
win_quote_trust_plot ## Figure B2

# Figure B3
win_trust <- lm(fav ~ win_negative * trust + win_qanon * trust, mturk) # favorability win, trust
win_trust_treat <- lm(fav ~ win_qanon * trust, mturk %>% filter(win_neutral != 1)) # favorability win, trust treatment comparison
win_trust_plot <- trust_me_plots(win_trust, win_trust_treat, 'win_negative', 'win_qanon', 1:4, 'Trust in Media') # plots
win_trust_plot ## Figure B3

# figure b4
# metro interaction
int_mod <- lm(therm ~ qanon * metro_area + wave, wc %>% filter(neutral == 0))
metro_me <- slopes(int_mod, variables = 'qanon', newdata = datagrid(metro_area = c(0,1)), conf_level = 1-(.05/2)) %>% 
  mutate(val = if_else(metro_area == 1, 'Metro', 'Non-Metro'))

metro_plot <- ggplot(metro_me, aes(x = val, y = estimate)) +
  geom_point() + 
  geom_errorbar(aes(x = val, ymin = conf.low, ymax = conf.high), width = 0) + 
  geom_hline(yintercept = 0, linetype = 'dashed') + 
  theme_baselike() + 
  theme(panel.border = element_rect(colour = "black", fill=NA, linewidth=1)) +
  labs(x = '', y = 'Marginal Effect of QAnon Support (vs Negative)') +
  coord_cartesian(ylim = c(-20,5))

# --------------------------------------------------------------------------- #

# Figure B5
# correlation between trust in fox and trust in media
cor(wc$trust_fox, wc$trust, use = 'pairwise.complete.obs')
# trust regression using trust in fox
trust_fox_mod <- lm(therm ~ negative * trust_fox + qanon * trust_fox, wc)
# marginal effects of trust in fox news
trust_fox_mes <- slopes(trust_fox_mod, variables = c('negative','qanon'), newdata = datagrid(trust_fox = 1:4), conf_level = 1-(.05/12)) %>% 
  mutate(comp = if_else(term == 'qanon', 'Q - Neu', 'Neg - Neu'))
# negative vs q comparison
trust_fox_comp <- lm(therm ~ qanon * trust_fox, wc %>% filter(neutral == 0))
trust_fox_comp_mes <- slopes(trust_fox_comp, variables = 'qanon', newdata = datagrid(trust_fox = 1:4), conf_level = 1-(.05/12)) %>% 
  mutate(comp = 'Q - Neg')

fox_df <- bind_rows(trust_fox_mes, trust_fox_comp_mes)
fox_df$comp <- factor(fox_df$comp, levels = c("Neg - Neu", "Q - Neu", "Q - Neg"))

(fox_trust_plot <- ggplot(fox_df, aes(x = as.factor(trust_fox), y = estimate)) + 
  geom_point(aes(shape = comp, color = comp), position = position_dodge(w = .3)) + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color = comp), width = 0, position = position_dodge(w = .3)) + 
  theme_baselike() + 
  geom_hline(yintercept = 0, linetype = 'dashed') + 
  labs(x = 'Trust Fox News', y = 'Marginal Effect of Treatment QAnon Support', color = 'Comparison', shape = 'Comparison') +
  scale_x_discrete(breaks = c(1:4), labels=c('None \nAt All', 'Not Very \nMuch', 'A Fair \nAmount', 'A Great \nDeal')) +
  scale_color_grey(name = 'Comparison', start = 0, end = .6) +
  scale_shape_manual(name = "Comparison", values = c(15, 16, 17)))

# --------------------------------------------------------------------------- #

# Figure B6, age and religion as alternative moderators 
age_mod <- lm(therm ~ qanon * age4 + wave, wc %>% filter(neutral == 0))
age_me <- slopes(age_mod, variables = 'qanon', newdata = datagrid(age4 = 1:4), conf_level = 1 - (.05/4))

evangelical_mod <- lm(therm ~ qanon * evangelical, wc %>% filter(neutral == 0))
ev_me <- slopes(evangelical_mod, variables = 'qanon', newdata = datagrid(evangelical = c(0,1)), conf_level = 1 - (.05/2))

rattend_mod <- lm(therm ~ qanon * relig_attend, wc %>% filter(neutral == 0))
relig_me <- slopes(rattend_mod, variables = 'qanon', newdata = datagrid(relig_attend = c(1,4,7,9)), conf_level = 1 - (.05/9))

alt_mod_mes <- bind_rows(
  age_me %>% mutate(label = 'Age', val = c('18-29', '30-44', '45-59', '60+')),
  ev_me  %>% mutate(label = 'Evangelical', val = c('Not Evangelical', 'Evangelical')),
  relig_me %>% mutate(label = 'Religious Attendance', val = c('Never', 'Several\nTimes\na Year', 'Nearly\nEvery\nWeek', 'Several\nTimes\na Week'))) %>% 
mutate(val = factor(val, 
  levels = c('18-29', '30-44', '45-59', '60+', 'Never', 'Not Evangelical', 'Evangelical', 'Several\nTimes\na Year', 'Nearly\nEvery\nWeek', 'Several\nTimes\na Week')))

alt_mods_plot <- ggplot(alt_mod_mes, aes(x = val, y = estimate)) +
  facet_wrap(~label, scales = 'free_x') + 
  geom_point() + 
  geom_errorbar(aes(x = val, ymin = conf.low, ymax = conf.high), width = 0) + 
  geom_hline(yintercept = 0, linetype = 'dashed') + 
  theme_baselike() + 
  theme(panel.border = element_rect(colour = "black", fill=NA, linewidth=1)) +
  labs(x = '', y = 'Marginal Effect of QAnon Support (vs Negative)') + 
  coord_cartesian(ylim = c(-50, 5))

# --------------------------------------------------------------------------- #

# Figure B7
# knowledge of qanon
know_q_mod <- lm(therm ~ qanon * dk_q + wave, wc %>% filter(neutral == 0))
know_q_mes <- slopes(know_q_mod, variables = 'qanon', newdata = datagrid(dk_q = c(0,1)), conf_level = 1 - (.05/2))

# importance of following news
newsimport_mod <- lm(therm ~ qanon * newsimport + wave, wc %>% filter(neutral == 0))
newsimport_mes <- slopes(newsimport_mod, variables = 'qanon', newdata = datagrid(newsimport = 1:5), conf_level = 1 - (.05/5))

both_know_mes <- bind_rows(know_q_mes %>% 
    mutate(label = 'Unfamilar with QAnon', 
      val = c('Support/Oppose QAnon','Unfamilar with QAnon')), 
  newsimport_mes %>% 
    mutate(label = 'Importance of Following News', 
      val = c(
        'Not at\nall important',
        'Not very\nimportant',
        'Somewhat\nimportant',
        'Very\nimportant',
        'Extremely\nimportant'))) %>% 
  mutate(val = factor(val, levels = c('Not at\nall important',
        'Not very\nimportant',
        'Somewhat\nimportant',
        'Very\nimportant',
        'Extremely\nimportant', 
        'Support/Oppose QAnon',
        'Unfamilar with QAnon')))

know_me_plot <- ggplot(both_know_mes, aes(x = val, y = estimate)) +
  facet_wrap(~label, scales = 'free_x') + 
  geom_point() + 
  geom_errorbar(aes(x = val, ymin = conf.low, ymax = conf.high), width = 0) + 
  geom_hline(yintercept = 0, linetype = 'dashed') + 
  theme_baselike() + 
  theme(panel.border = element_rect(colour = "black", fill=NA, linewidth=1)) +
  labs(x = '', y = 'Marginal Effect of QAnon Support (vs Negative)')

# additional robustness checks
# Figure B9
wc <- wc %>% 
  mutate(trust2 = case_when(trust %in% 1:2 ~ 0, trust %in% 3:4 ~ 1),
    gop = if_else(dem == 0 & ind == 0, 1, 0))

# trust binary
ols_trust2 <- lm(therm ~ negative * trust2 + qanon * trust2 + wave, wc) 
ols_trust2_treat <- lm(therm ~ qanon * trust2 + wave, wc %>% filter(neutral != 1)) # therm ~ qanon-negative x trust

trust2_me1 <- slopes(ols_trust2, variables = c('negative', 'qanon'), newdata = datagrid(trust2 = c(0,1)), conf_level = 1 - (.05/6)) %>% 
  mutate(term2 = c(rep('Neg – Neu', 2), rep('Q — Neu', 2)))
trust2_me2 <- slopes(ols_trust2_treat, variables = c('qanon'), newdata = datagrid(trust2 = c(0,1)), conf_level = 1 - (.05/6)) %>% 
  mutate(term2 = 'Q – Neg')

trust2_df <- bind_rows(trust2_me1, trust2_me2) %>% 
  mutate(trust2 = if_else(trust2 == 1, 'High Trust in Media', 'Low Trust in Media'))

trust2_plot <- ggplot(trust2_df, aes(x = trust2, y = estimate)) + 
  geom_point(aes(shape = term2, color = term2), position = position_dodge(w = .3)) + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color = term2), width = 0, position = position_dodge(w = .3)) + 
  theme_baselike() + 
  geom_hline(yintercept = 0, linetype = 'dashed') + 
  labs(x = 'Trust in Media', y = 'Marginal Effect of Treatment', color = 'Comparison', shape = 'Comparison') +
  scale_color_grey(name = 'Comparison', start = 0, end = .6) +
  scale_shape_manual(name = "Comparison", values = c(15, 16, 17))

# pid binary
ols_gop <- lm(therm ~ negative * gop + qanon * gop + wave, wc) 
ols_gop_treat <- lm(therm ~ qanon * gop + wave, wc %>% filter(neutral != 1)) # therm ~ qanon-negative x trust

gop_me1 <- slopes(ols_gop, variables = c('negative', 'qanon'), newdata = datagrid(gop = c(0,1)), conf_level = 1 - (.05/6)) %>% 
  mutate(term2 = c(rep('Neg – Neu', 2), rep('Q — Neu', 2)))
gop_me2 <- slopes(ols_gop_treat, variables = c('qanon'), newdata = datagrid(gop = c(0,1)), conf_level = 1 - (.05/6)) %>% 
  mutate(term2 = 'Q – Neg')

gop_df <- bind_rows(gop_me1, gop_me2) %>% 
  mutate(gop2 = if_else(gop == 1, 'Republican', 'Democrat/Independent'))

gop_plot <- ggplot(gop_df, aes(x = gop2, y = estimate)) + 
  geom_point(aes(shape = term2, color = term2), position = position_dodge(w = .3)) + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color = term2), width = 0, position = position_dodge(w = .3)) + 
  theme_baselike() + 
  geom_hline(yintercept = 0, linetype = 'dashed') + 
  labs(x = 'Party', y = 'Marginal Effect of Treatment', color = 'Comparison', shape = 'Comparison') +
  scale_color_grey(name = 'Comparison', start = 0, end = .6) +
  scale_shape_manual(name = "Comparison", values = c(15, 16, 17))

# make figure B9
bin_var_plot <- ggarrange(trust2_plot, gop_plot, ncol = 1)

# Figure B10
ols_ideo <- lm(therm ~ negative * ideology + qanon * ideology + wave, wc) 
ols_ideo_treat <- lm(therm ~ qanon * ideology + wave, wc %>% filter(neutral != 1)) 

ideo_mes <- bind_rows(
  slopes(ols_ideo, variables = c('negative', 'qanon'), newdata = datagrid(ideology = c(1:5)), conf_level = 1 - (0.05/15)) %>% 
    mutate(term2 = c(rep('Neg — Neu', 5), rep('Q — Neu', 5))),
  slopes(ols_ideo_treat, variables = c('qanon'), newdata = datagrid(ideology = c(1:5)), conf_level = 1 - (0.05/15)) %>% 
    mutate(term2 = 'Q — Neg')
  ) %>% mutate(term2 = factor(term2, levels = c('Neg — Neu', 'Q — Neu', 'Q — Neg')))

ideo_plot <- ggplot(ideo_mes, aes(x = ideology, y = estimate)) + 
  geom_point(aes(shape = term2, color = term2), position = position_dodge(w = .3)) + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color = term2), width = 0, position = position_dodge(w = .3)) + 
  theme_baselike() + 
  geom_hline(yintercept = 0, linetype = 'dashed') + 
  labs(x = 'Ideology (Liberal to Conservative)', y = 'Marginal Effect of Treatment', color = 'Comparison', shape = 'Comparison') +
  scale_color_grey(name = 'Comparison', start = 0, end = .6) +
  scale_shape_manual(name = "Comparison", values = c(15, 16, 17))